module mymodule

use select_precision
use constants
use complete_elliptic
implicit none

contains

    pure function CircularCoaxImpedance( d_outer , d_inner , t , epsilon_r )
        !Calculates the impedance of an offset cicular offset coaxial line from 
        ! the analytical solution
        implicit none
        real(kind=wp),intent(in)    :: d_outer, d_inner , t , epsilon_r
        real(kind=wp)                   :: CircularCoaxImpedance , phi  , acosh
        
        if ( t+HALF*d_inner >= HALF*d_outer ) then
            CircularCoaxImpedance = ZERO
            return
        elseif   ( t < epsilon(t) ) then
            CircularCoaxImpedance = HALF*PIINV*ETA0/sqrt(epsilon_r)*log( d_outer/d_inner )
            return
        else        
            phi = TWO*t/d_outer
            CircularCoaxImpedance = HALF*PIINV*ETA0/sqrt(epsilon_r)*&
                acosh(  HALF*d_outer/d_inner*(ONE-phi*phi) + HALF*d_inner/d_outer )
            return
        endif
    end function CircularCoaxImpedance
    
    pure function Stripline( b , w , epsilon_r )
        implicit none
        real(kind=wp),intent(in)    :: b,w ,  epsilon_r
        real(kind=wp)                   :: Stripline , k_squared , k_p_squared
        
        k_squared = ONE/cosh( PIHALF*w/b )**2
        k_p_squared =  ONE - k_squared
        Stripline = QUARTER*ETA0/sqrt(epsilon_r)*CompleteEllipticFirstKind( k_p_squared )/CompleteEllipticFirstKind( k_squared )
        return
    end function Stripline
    
    pure function linspace( xstart , xstop , N )
        implicit none
        real(kind=wp),intent(in)    :: xstart,xstop
        integer,intent(in)          :: N
        integer                     :: i

        real(kind=wp),allocatable,dimension(:)  :: linspace
        
        allocate( linspace(N) )
        
        
        do i=1,N
            linspace(i) = (xstop-xstart) * real(i-1,wp) / real(N-1,wp) + xstart
        enddo
        
        return
    end function linspace
    
    pure function SquareCoaxImpedance( d_outer , d_inner , epsilon_r )
        !Calculates the impedance of a square coaxial line given x and epsilon_r
        ! x is the solution to the lambda relation which must be found prior to calling this function
        implicit none
        real(kind=wp),intent(in) :: d_outer , d_inner , epsilon_r
        real(kind=wp)                :: SquareCoaxImpedance  , k_prime , k , ratio , m  , m_prime
        
        ratio = (d_outer-d_inner)/(d_outer+d_inner)
        
        m = find_zero( ZERO , ONE , ONE*epsilon(ZERO) , lambda_relation , ratio )
        !m is now the value such that K( x' )/K( x ) = ratio ; x'=sqrt(1-x**2)
        
        m_prime = ONE - m
        m = (sqrt(m_prime) - sqrt(m))**4/(sqrt(m_prime) + sqrt(m))**4
        m_prime =  ONE - m

        SquareCoaxImpedance = ONEEIGTH*ETA0/sqrt(epsilon_r)* &
            CompleteEllipticFirstKind( m )/CompleteEllipticFirstKind( m_prime )
            
        return
    end function SquareCoaxImpedance

    pure function lambda_relation( m ,ratio) 
        implicit none
        real(kind=wp),intent(in) :: m,ratio
        real(kind=wp)            :: lambda_relation , m_prime
        ! This function is a bit convolved but basically we seek the value x such that
        ! K( x' )/K( x ) = ratio
        ! where x' = sqrt( 1-x**2 ) and ratio is the ratio as defined by the square coaxial geometry
        ! Remembering that the elk function takes the complementary m_c = 1-m instead of k,
        ! that m = k**2 and that m_prime = 1-m we arrive at the solution below
        
        m_prime = ONE - m ! = k_prime**2
        
        lambda_relation = CompleteEllipticFirstKind( m )/CompleteEllipticFirstKind( m_prime ) - ratio
        return
    end function lambda_relation

    pure function find_zero ( a, b, t, func, second_argument )
!*****************************************************************************80
!
!! ZERO seeks the root of a function Func(X) in an interval [A,B].
!
!  Discussion:
!
!    The interval [A,B] must be a change of sign interval for F.
!    That is, F(A) and F(B) must be of opposite signs.  Then
!    assuming that F is continuous implies the existence of at least
!    one value C between A and B for which F(C) = 0.
!
!    The location of the zero is determined to within an accuracy
!    of 6 * MACHEPS * abs ( C ) + 2 * T.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    12 April 2008
!
!  Author:
!
!    Original FORTRAN77 version by Richard Brent.
!    FORTRAN90 version by John Burkardt.
!    Fortran 95 version by Magnus Lundmark
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization Without Derivatives,
!    Dover, 2002,
!    ISBN: 0-486-41998-3,
!    LC: QA402.5.B74.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, B, the endpoints of the change of 
!    sign interval.
!
!
!    Input, real ( kind = 8 ) T, a positive error tolerance.
!
!    Input, external real ( kind = 8 ) F, the name of a user-supplied
!    function, of the form "FUNCTION F ( X )", which evaluates the
!    function whose zero is being sought.
!
!    Output, real ( kind = 8 ) ZERO, the estimated value of a zero of
!    the function F.
!
    implicit none
    real(kind = wp), intent(in) :: a,b,t,second_argument
    real ( kind = wp )                 :: c,d,e,fa,fb,fc,m,machep,p,q,r,s,sa,sb,tol
    real ( kind = wp )                 ::find_zero
    interface
        pure function func(x,y)
            use select_precision
            implicit none
            real(wp), intent(in)  :: x,y
            real(wp)              :: func
        end function func
    end interface
  
    machep = epsilon( a )
!
!  Make local copies of A and B.
!
    sa = a
    sb = b
    fa = func ( sa ,second_argument)
    fb = func ( sb ,second_argument)

    c = sa
    fc = fa
    e = sb - sa
    d = e

    do

        if ( abs ( fc ) < abs ( fb ) ) then
            sa = sb
            sb = c
            c = sa
            fa = fb
            fb = fc
            fc = fa
        endif

        tol = TWO * machep * abs ( sb ) + t
        m = HALF * ( c - sb )

        if ( abs ( m ) <= tol .or. fb == ZERO ) then
            exit
        endif

        if ( abs ( e ) < tol .or. abs ( fa ) <= abs ( fb ) ) then
            e = m
            d = e
        else
            s = fb / fa
            if ( sa == c ) then
                p = TWO * m * s
                q = ONE - s
            else
                q = fa / fc
                r = fb / fc
                p = s * ( TWO * m * a * ( q - r ) - ( sb - sa ) * ( r - ONE ) )
                q = ( q - ONE ) * ( r - ONE ) * ( s - ONE )
            endif
            
            if ( ZERO < p ) then
                q = - q
            else
                p = - p
            endif
            s = e
            e = d

            if ( TWO * p < THREE * m * q - abs ( tol * q ) .and. p < abs ( HALF * s * q ) ) then
                d = p / q
            else
                e = m
                d = e
            endif

        endif

        sa = sb
        fa = fb

        if ( tol < abs ( d ) ) then
            sb = sb + d
        else if ( ZERO < m ) then
            sb = sb + tol
        else
            sb = sb - tol
        endif

        fb = func( sb ,second_argument )

        if ( ( ZERO < fb .and. ZERO < fc ) .or. ( fb <= ZERO .and. fc <= ZERO ) ) then
            c = sa
            fc = fa
            e = sb - sa
            d = e
        endif

    enddo

    find_zero = sb
    return
    end function find_zero

end module mymodule